Northeastern US Shelf Region 2023: Sea Surface Temperature
This report was created to track the sea surface temperature regimes for the Northeastern US continental shelf, an area sampled as part of the NMFS Northeast trawl survey.
DISCLAIMER: Any data within 2-weeks of the current date are subject to revision and may change. Please use caution when reporting information that contains these values.
Northeastern US Shelf Region On a Map
Whenever discussing the Northeastern US Shelf Region in this report, we refer to the following spatial extent displayed below. The coordinates of this bounding box are the same coordinates used to clip the sea surface temperature data.
1: Temperature History
Using the coordinates shown above we can create an area-specific temperature history. For any day of the year since September of 1981, data is available to calculate the average sea surface temperature within that area for each day.
Additionally, we can compare observed temperatures against the expected conditions based on a climatology using a specified reference period. The standard reference period used for the climatology here is 1982-2011, a 30-year period.
Climatological Patterns of the Northeastern US Shelf Region
Daily climatologies are a record of daily averages based on the day of the year. They record what the average temperature has been for each day of the year, across the specified reference period.
In addition to the daily average, we also look at how variable temperatures are. This variability is used to benchmark how “rare” extreme events are to determine whether they are part of a natural cycle or not.
Common benchmarks for extreme events are the 10th and 90th percentile. Temperatures above the 90th percentile are among the hottest 10% of days in the reference period. Temperatures below the 10th percentile correspond to the coldest 10%. The remaining 80% of days fall somewhere in-between and showcase the range of temperatures we might expect to occur given the natural variability of the climate.
Working in Anomalies
Using the climatology as a reference, we can see hot much “hotter” or “colder” a given day is than what we would on average expect. This difference from the average is what we call anomalies.
In the table below, Sea Surface Temperature is the mean temperature observed for that date averaged across all cells within the area. Climate Avg. & Climate SD are the climate means and standard deviations for a 1982-2011 climatology. Temperature Anomaly is the daily observed sea surface temperature minus the climate mean.
In this way a Temperature Anomaly is just: How much hotter or colder is it than the average temperature for that day.
Northeastern US Shelf Region - Regional Sea Surface Temperature
Temperature Unit: Celsius
Date
Sea Surface Temperature
Day of Year
Climate Avg.
Temperature Anomaly
2022-12-26
10.13
361
8.77
1.36
2022-12-27
9.88
362
8.73
1.14
2022-12-28
9.38
363
8.68
0.70
2022-12-29
9.10
364
8.60
0.50
2022-12-30
9.02
365
8.51
0.51
2022-12-31
9.09
366
8.41
0.67
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Climatology Reference Period: 1982-2011.
Detecting Long-Term Changes
One reason we go through this trouble of calculating climatologies and anomalies is to set expectations around how much we expect things to naturally vary, and to detect when things fall outside that range.
In the case of the Northeastern US Shelf Region we are now consistently outside of the range of temperatures normally expected from 1982-2011. The following plot colors the monthly average by how far temperatures are from the climatological average. Blue colors represent a month that was cooler than average, with red indicating a warmer than average month.
Warming Rates
Seeing this change towards a warmer climate, it is then natural to ask how rapidly the change is happening. This is where we turn to warming rates. The warming trends below were calculated using all the available data for complete years beginning with 1982 through the end of 2022.
The overlaid trend lines then track how warming has increased with time. A dotted line has been included to show how the global average temperature has changed during the same period.
Warming Trends
Shift in Regional Ocean Circulation
Beginning around 2008, temperature in the region swung upward, and since 2010 most yearly temperatures have been over 1 degree above the hematological average of the previous 30-years. With some years experiencing temperatures more than 2 degrees C above that average.
The magnitude of this regime shift was 1.09C, or 1.96F. Bringing average yearly temperatures up from the 30-year average of 11.87C to 12.96C, or 53.37F to 55.33F.
Seasonal Patterns
Overall Temperature Increase
Marine Heatwaves
A marine heatwave is defined as a situation when seawater temperatures exceeds a seasonally-varying threshold (usually the 90th percentile) for at least 5 consecutive days. Successive heatwaves with gaps of 2 days or less are considered part of the same event. The heatwave threshold used below was 90%. The heatwave history for Northeastern US Shelf Region is displayed below:
Note: For the figures below heatwave events were determined using the methods of Hobday et al. 2016 and implemented using the R package {heatwaveR}.
Heatwave Events
Heatwave Trends
2: 2023 Observations
Temperatures of the Last Full Year (Interactive)
Daily temperatures for the last 365 days are displayed below with reference to how they fall against the marine heatwave and cold spell thresholds.
2023 Temperatures
Temperatures for the current year can be seen against the same thresholds with the following plot:
2023 Anomalies
Changing from absolute temperatures to anomalies reveals the degree to which this year is departure from what we would have expected from the climatology of 1982-2011’s temperatures.
2023 in Context
When 2023 is overlaid against the daily values and climatological average for the region, we can see how temperatures compare to the history of the region.
3: Global Temperature Maps
The 2022 global sea surface temperature anomalies have been loaded and displayed below to visualize how different areas of the ocean experience swings in temperature.
For perspective on where excess heat is distributed around the world, here first are maps of anomalies at a global perspective:
Annual Average
Winter
Spring
Summer
Fall
4: Regional Temperature Maps
The regional views of 2022’s sea surface temperature anomalies have been loaded and displayed below to visualize how localized differences changed throughout the year.
Annual Average
Winter
Spring
Summer
Fall
Heatwave Progression
Looking specifically at the last heatwave event, we can step through how the event progressed over time, and developing pockets or warmer/colder water masses.
5: Ranking Warming Rates
If we look at the rates of change from 1982-2022 for each grid cell, rather than the observed temperature, it is possible to rank how hot each location on earth is warming relative to all the others.
Once we have the rankings, we can then take the average ranking within the Gulf of Maine we can obtain the average warming rank for the area compared to the rest of the globe.
Based on data from 1982-2022, the warming rates of Northeastern US Shelf Region have been some of the highest in the world. The area as a whole has been increasing at a rate of 0.04\(^{\circ}C/year\) which is faster than 92.1% of the world’s oceans.
Over that same period locations within the Northeastern US Shelf Region have been warming at rates as low as -0.016\(^{\circ}C/year\) and as rapidly as 0.074\(^{\circ}C/year\), corresponding to ranks as low as 0.3% and as high as 99.9%.
Mapped below are the corresponding warming rates and their global rankings.
Warming Rates
Warming Ranks
6: Shifting Baselines
In 2022 NOAA is transitioning standard climatologies from the 30-year period of 1982-2011 to a new period spanning 1992-2020. Changes in climate regimes often does not result in a uniform upward or downward change that is consistent throughout the year.
The plot below shows just how both the average temperature, as well as the annual highs and lows have shifted. When looking specifically at Northeastern US Shelf Region here is how the expected temperature for each day of the year has shifted.
From this we can see that the Fall temperatures have risen more than those of the spring. There is also a large change in where the threshold for Marine Heatwave events sits, a consequence of exceptionally warm Fall temperatures becoming more common.
---title: "Northeastern US Continental Shelf SST Report"author: name: "Adam Kemberling" url: https://github.com/adamkemberling affiliation: Gulf of Maine Research Institutedescription: | Temperature tracking for the Northeastern US Continental Shelf Regiondate: "Updated on: `r Sys.Date()`"format: html: code-fold: true code-tools: true df-print: kable self-contained: trueeditor: visualexecute: echo: false warning: false message: false fig.height: 6 fig.width: 8 fig.align: "center" comment: ""---```{r}# Packageslibrary(lubridate) # date supportlibrary(raster) # raster supportlibrary(rnaturalearth) # coastline polygonslibrary(sf) # simple feature supportlibrary(stars) # raster plotting with ggplotlibrary(gmRi) # styling for GMRIlibrary(scales) # axis labelslibrary(here) # project navigation# library(janitor) # data cleaninglibrary(patchwork) # multiplot arrangementlibrary(tidyverse) # data wrangling and plottinglibrary(heatwaveR) # heatwave detectionlibrary(plotly) # interactive plotslibrary(ggpmisc) # statistics added to plotslibrary(gt) # table stylinglibrary(ggalt) # Geom dumbbell and lollipoplibrary(geomtextpath) # text within geom_lines or geom_smooth library(xaringanExtra) # tab panelslibrary(ggforce) # smart text labels in plots & facet_zoomlibrary(ggthemes) # theme supportlibrary(ggrepel) # text labelinglibrary(ggtext) # Colored ggplot textlibrary(xaringanExtra) # For tabsets# Support FunctionssuppressPackageStartupMessages(source(here("R/oisst_support_funs.R")))suppressPackageStartupMessages(source(here("R/temp_report_support.R")))# File Pathsmills_path <-cs_path("mills")oisst_path <-cs_path("res", "OISST/oisst_mainstays")# Set ggplot theme for figurestheme_set(theme_bw() +theme(# Titlesplot.title =element_text(hjust =0, face ="bold", size =14),plot.subtitle =element_text(size =9),plot.caption =element_text(size =7.2, margin =margin(t =20), color ="gray40"),legend.title =element_text(size =9),legend.text =element_text(size =7.5),# Axesaxis.line.y =element_line(),axis.ticks.y =element_line(), axis.line.x =element_line(),axis.ticks.x =element_line(), axis.text =element_text(size =11),# Facetsstrip.text =element_text(color ="white", face ="bold",size =11),strip.background =element_rect(color ="white", fill ="#00736D", linewidth =1, linetype="solid") ) )# Polygons for mappingnew_england <-ne_states("united states of america", returnclass ="sf")canada <-ne_states("canada", returnclass ="sf")world <-ne_countries(returnclass ="sf")greenland <-ne_states(country ="greenland", returnclass ="sf")# Adding Logo to plotlogo_path <-paste0(system.file("stylesheets", package ="gmRi"), "/gmri_logo.png")lab_logo <- magick::image_read(logo_path)``````{r params formatting}# Set the shapefile to use as a parameterparams <-list(region ="inuse_strata")# File paths for various extents based on params$regionregion_paths <-get_timeseries_paths(region_group ="nmfs_trawl_regions", box_location ="cloudstorage")# clean up the nametidy_name <-"Northeastern US Shelf Region"# Current Yearcurrent_yr <-year(Sys.Date())# Use panelsetuse_panelset()# Turn on the stylesheet:gmRi::use_gmri_style_rmd(css ="gmri_rmarkdown.css")# Turn on the gmri font for plotsshowtext::showtext_auto()```# `r paste(tidy_name, current_yr)`: Sea Surface TemperatureThis report was created to track the sea surface temperature regimes for the Northeastern US continental shelf, an area sampled as part of the NMFS Northeast trawl survey.Satellite sea surface temperature data used was obtained from the National Center for Environmental Information (NCEI). With all maps and figures displaying [NOAA's Optimum Interpolation Sea Surface Temperature Data](https://www.ncdc.noaa.gov/oisst).**DISCLAIMER:** Any data within 2-weeks of the current date are subject to revision and may change. Please use caution when reporting information that contains these values.### `r tidy_name` On a MapWhenever discussing the `r tidy_name` in this report, we refer to the following spatial extent displayed below. The coordinates of this bounding box are the same coordinates used to clip the sea surface temperature data.```{r}# Load the bounding box for Northeast Shelfpoly_path <- region_paths[[params$region]][["shape_path"]]region_extent <-st_read(poly_path, quiet =TRUE)# region_extent <- st_union(region_extent)# Pull extents for the region to set crop extentcrop_x <-st_bbox(region_extent)[c(1,3)]crop_y <-st_bbox(region_extent)[c(2,4)]# Map the study areamap_study_area(region_extent =st_union(region_extent),x_buffer =c(-2, 2),y_buffer =c(-0.75, 0.75), new_england_sf = new_england, canada_sf = canada,greenland_sf = greenland,shape_linetype =2)```# 1: Temperature HistoryUsing the coordinates shown above we can create an area-specific temperature history. For any day of the year since **September of 1981**, data is available to calculate the average sea surface temperature within that area for each day.Additionally, we can compare observed temperatures against the expected conditions based on a climatology using a specified reference period. The standard reference period used for the climatology here is 1982-2011, a 30-year period.```{r load data}# Load the timeseriestimeseries_path <- region_paths[[params$region]][["timeseries_path"]]region_timeseries <-read_csv( timeseries_path, col_types =cols(), guess_max =1e6) %>%mutate(tod =format(time, format ="%H:%M:%S")) # Clean up the data - add labelsregion_timeseries <- region_timeseries %>%mutate(time =as.Date(time)) %>%distinct(time, .keep_all = T) %>%supplement_season_info() %>%filter(year %in%c(1982:2022))# Summarize by year to return mean annual anomalies and varianceannual_summary <- region_timeseries %>%filter(year %in%c(1982:2022)) %>%group_by(year) %>%temperature_summaries() %>%mutate(yr_as_dtime =as.Date(paste0(year, "-07-02")),anom_direction =ifelse(area_wtd_anom >0, "Hot", "Cold"))# Get heatwave statuses for each day:# Uses area weighted sst by defaultregion_hw <-pull_heatwave_events(temperature_timeseries = region_timeseries, threshold =90, clim_ref_period =c("1982-01-01", "2011-12-31")) %>%supplement_hw_data() %>%filter(doy !=366)# Last Calendar Yearthis_yr <- region_hw %>%filter(year(time) ==2022)# Set up the date within a year for Heatmapbase_date <-as.Date("2000-01-01")grid_data <- region_hw %>%mutate(year =year(time),yday =yday(time),flat_date =as.Date(yday-1, origin = base_date),yr_rev =factor(year),yr_rev =fct_rev(yr_rev))# # Global Temperature Anomaly Ratesglobal_anoms <-read_csv(paste0(oisst_path, "global_timeseries/global_anoms_1982to2011.csv"), guess_max =1e6,col_types =cols()) %>%mutate(year =year(time)) %>%filter(between(year, 1982, 2022))# summarize by year againglobal_summary <- global_anoms %>%group_by(year) %>%temperature_summaries() %>%mutate(yr_as_dtime =as.Date(paste0(year, "-07-02")),anom_direction =ifelse(area_wtd_anom >0, "Hot", "Cold"))```## Climatological Patterns of the `r tidy_name`Daily climatologies are a record of daily averages based on the day of the year. They record what the average temperature has been for each day of the year, across the specified reference period.In addition to the daily average, we also look at how variable temperatures are. This variability is used to benchmark how "rare" extreme events are to determine whether they are part of a natural cycle or not.Common benchmarks for extreme events are the 10th and 90th percentile. Temperatures above the 90th percentile are among the hottest 10% of days in the reference period. Temperatures below the 10th percentile correspond to the coldest 10%. The remaining 80% of days fall somewhere in-between and showcase the range of temperatures we might expect to occur given the natural variability of the climate.```{r}clim_colors <-c("Daily Temperature"="gray80","Climatological Average"="gray10","90th Percentile"="coral3","10th Percentile"="royalblue")# pull one year for the lines, the reference period for pointsone_yr <- grid_data %>%filter(yr ==2022)ref_yrs <- grid_data %>%filter(yr %in%c(1982:2011))# Make a plotggplot(data = one_yr, aes(x = flat_date)) +geom_point(data = ref_yrs, aes(y = sst, color ="Daily Temperature"), alpha =0.6) +geom_line(aes(y = seas, color ="Climatological Average"), size =1) +geom_textline(aes(y = mhw_thresh, color ="90th Percentile", label ="Heatwave Threshold"), linetype =1, size =4, hjust =0.4, linewidth =1) +geom_textline(aes(y = mcs_thresh, color ="10th Percentile", label ="Coldspell Threshold"), linetype =1, size =4, hjust =0.6, linewidth =1) +scale_x_date(date_labels ="%b", date_breaks ="1 month", expand =expansion(add =c(0,0))) +scale_y_continuous(labels =number_format(suffix =" \u00b0C")) +scale_color_manual(values = clim_colors) +labs(x ="Day of the Year", y ="Sea Surface Temperature",title =str_c(tidy_name, " Daily Climatology"),subtitle ="Anticipated Daily Temperatures Based on a 1982-2011 Climate") +#theme_gmri() +theme(legend.position ="bottom",legend.key.height =unit(.5, "lines"),legend.key.width =unit(5, "lines"),plot.margin =unit(c(.5,1,.3,.5), "cm"),plot.title.position ="plot",legend.margin =margin(c(7,0,0,0)),legend.justification ="center") +guides(color =guide_legend(title ="",title.hjust =0.5,nrow =1,title.position ="left",label.position ="top", override.aes =list(color =c("gray50", "gray10", "coral3", "royalblue"),shape =c(16, NA, NA, NA), linetype =c(0, 1, 2, 2),alpha =c(1,1,1,1), label =c("", "", "", "")),byrow = T))```### Working in AnomaliesUsing the climatology as a reference, we can see hot much "hotter" or "colder" a given day is than what we would on average expect. This difference from the average is what we call anomalies.In the table below, `Sea Surface Temperature` is the mean temperature observed for that date averaged across all cells within the area. `Climate Avg.` & `Climate SD` are the climate means and standard deviations for a 1982-2011 climatology. `Temperature Anomaly` is the daily observed sea surface temperature minus the climate mean.In this way a `Temperature Anomaly` is just: How much hotter or colder is it than the average temperature for that day.```{r}# Display Table of last 6 entriestail(region_timeseries) %>%mutate(across(where(is.numeric), round, 2)) %>%select(Date = time,`Sea Surface Temperature`= sst,#`Area-Weighted SST` = area_wtd_sst,`Day of Year`= modified_ordinal_day,`Climate Avg.`= sst_clim,#`Area-weighted Climate` = area_wtd_clim,`Temperature Anomaly`= sst_anom#,#`Area-Weighted Anomaly` = area_wtd_anom ) %>%gt() %>%tab_header(title =md(paste0("**", tidy_name, " - Regional Sea Surface Temperature", "**")), subtitle =paste("Temperature Unit: Celsius")) %>%tab_source_note(source_note =md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>%tab_source_note(md("*Climatology Reference Period: 1982-2011.*"))# march 1st sstmar1 <- region_timeseries %>%filter(modified_ordinal_day ==61) %>%distinct(sst_clim) %>%pull(sst_clim)```### Detecting Long-Term ChangesOne reason we go through this trouble of calculating climatologies and anomalies is to set expectations around how much we expect things to naturally vary, and to detect when things fall outside that range.In the case of the `r tidy_name` we are now consistently outside of the range of temperatures normally expected from 1982-2011. The following plot colors the monthly average by how far temperatures are from the climatological average. Blue colors represent a month that was cooler than average, with red indicating a warmer than average month.```{r, temp stripes, fig.height=2.25}# Make Monthly Averagesmonthly_temps <- region_timeseries %>%group_by(year =year(time),month =month(time)) %>%summarise(sst =mean(area_wtd_sst),anom =mean(area_wtd_anom),.groups ="drop") %>%mutate(time =as.Date(str_c(year, str_pad(month, width =2, pad ="0", side ="left"), "01", sep ="-")))# Plot monthly timelineggplot(monthly_temps, aes(x = time, y =0)) +geom_tile(aes(fill = anom)) +scale_fill_distiller(palette ="RdBu", limits =c(-1.5,1.5), oob = oob_squish) +scale_y_continuous(expand =expansion(add =c(0,0.01))) +scale_x_date(expand =expansion(add =c(45,45))) +theme(legend.position ="none",plot.title =element_text(hjust =0.5, face ="bold", size =14),plot.subtitle =element_text(hjust =0.125),panel.border =element_rect(color ="black", size =1, fill ="transparent"),axis.ticks.y =element_blank(),axis.text.y =element_blank(),axis.title.y =element_blank(),panel.grid =element_blank(),plot.caption =element_text(size =7.2, margin =margin(t =20), color ="gray40") ) +labs(title =str_c(tidy_name, "\nMonthly Sea Surface Temperature Anomalies"),x ="",caption ="Monthly Averages from 1982-2022 | Climatology Period 1982-2011")```### Warming RatesSeeing this change towards a warmer climate, it is then natural to ask how rapidly the change is happening. This is where we turn to warming rates. The warming trends below were calculated using all the available data for complete years beginning with 1982 through the end of 2022.The overlaid trend lines then track how warming has increased with time. A dotted line has been included to show how the global average temperature has changed during the same period.#### Warming Trends```{r annual regressions}#### Run Warming Rate regressions ####rate_data <-list("shelf F"=get_decadal_rates(temp_df = annual_summary, temp_col ="area_wtd_f", year_col ="year", year_lim =c(1982, 2022), area_name ="Northeast US Shelf", degree_c = F),"shelf C"=get_decadal_rates(temp_df = annual_summary, temp_col ="area_wtd_sst", year_col ="year", year_lim =c(1982, 2022), area_name ="Northeast US Shelf", degree_c = T),"Global All"=get_decadal_rates(temp_df = global_summary, temp_col ="area_wtd_f", year_col ="year", year_lim =c(1982, 2022), area_name ="Global", degree_c = F),"Global All C"=get_decadal_rates(temp_df = global_summary, temp_col ="area_wtd_sst", year_col ="year", year_lim =c(1982, 2022), area_name ="Global", degree_c = T))# Build Regression Equation Labelseq_all_c <- rate_data[["shelf C"]][["eq_label"]]eq_all <- rate_data[["shelf F"]][["eq_label"]]eq_15 <- rate_data[["shelf 15"]][["eq_label"]]eq_global <- rate_data[["Global All"]][["eq_label"]]eq_global_c <- rate_data[["Global All C"]][["eq_label"]]# Generate a smoothed temperature line using splinesyearly_temp_smooth <-as.data.frame(spline(annual_summary$yr_as_dtime, annual_summary$area_wtd_anom)) %>%mutate(x =as.Date(x, origin ="1970-01-01"))``````{r annual trend plot}# Fancy colored markdown titles and subtitles:neshelf_rate <-str_c(rate_data$`shelf C`$decadal_rate /10, "\u00b0C / year")glob_rate <-str_c(rate_data$`Global All C`$decadal_rate /10, "\u00b0C / year")fancy_title <-str_c("<span style='color:#00608A'>",tidy_name, "</span>", " Warming Faster than <span style='color:#407331'>Global Avg.</span>")fancy_subtitle <-str_c(tidy_name, "'s annual warming rate of <span style='color:#00608A'>",neshelf_rate, "</span>, is 2.5X the global average of <span style='color:#407331'>", glob_rate,"</span>.")#### Annual Trend Plot ####ggplot(data = annual_summary, aes(yr_as_dtime, area_wtd_anom)) +# Add daily datageom_line(data = region_timeseries,aes(time, area_wtd_anom, color ="Daily Temperatures")) +# Overlay yearly means# geom_line(data = yearly_temp_smooth, aes(x, y, color = "Average Yearly Temperature"), alpha = 0.8, linetype = 2) +geom_line(aes(color ="Average Yearly Temperature"), alpha =0.8, linetype =1, linewidth =1) +geom_point(color ="gray10", alpha =0.7, size =0.75) +# Add regression linesgeom_textsmooth(data =filter(global_summary, year <=2022),method ="lm", text_smoothing =30,label ="Global Trend",color =gmri_cols("green"),linewidth =1,formula = y ~ x, se = F,linetype =1, hjust =0.925) +geom_smooth(data =filter(annual_summary, year <=2022),method ="lm",aes(color ="1982-2022 NE Shelf Trend"), formula = y ~ x, se = F,linetype =1) +# Colorsscale_color_manual(values =c("1982-2022 NE Shelf Trend"=gmri_cols("gmri blue"),"Average Yearly Temperature"="gray10","Daily Temperatures"="gray90")) +# Axesscale_y_continuous(sec.axis =sec_axis(trans =~as_fahrenheit(., data_type ="anomalies"),labels =number_format(suffix =" \u00b0F")),labels =number_format(suffix =" \u00b0C")) +# labels + themelabs(x ="", y ="Sea Surface Temperature Anomaly",caption =paste0("Anomalies calculated using 1982-2011 reference period."),title = fancy_title,subtitle = fancy_subtitle) +# themetheme(plot.title =element_markdown(),plot.subtitle =element_markdown(),legend.title =element_blank(),legend.position ="bottom",legend.background =element_rect(fill ="transparent"),legend.key =element_rect(fill ="transparent", color ="transparent")) +guides(color =guide_legend(title ="",title.hjust =0.5,nrow =1,title.position ="left",label.position ="top", override.aes =list(point =c(0, 0, 0),linetype =c(1, 1, 1),byrow = T)))``````{r annual trend plot simplified, eval = FALSE}# Fahrenheit plottemp_simplified <-global_rate_comparison(annual_summary_dat = annual_summary, global_summary_dat = global_summary, eq_all =str_remove(eq_all, " 1982-2022"), eq_global =str_remove(eq_global, " 1982-2022"), temp_units ="F",region_label ="Northeast US Shelf") +coord_cartesian(clip ="off") +geom_logo(x_npc =0.115, y_npc =-0.135, logo_height =1, height_units ="cm")# add logotemp_simplified ``````{r}#| label: solid-lane-version# Get the global temperature warming rate# Decadal Rateglobal_yr_mod <-lm(anom_f ~ year, global_summary)global_annual_rate <-coef(global_yr_mod)[[2]]global_decadal_rate_f <-round((global_annual_rate*10), 2)# Join the data we need together,# Put the description into the dataframe so we can call it using one geom callglob_plot_dat <-bind_rows(list(# 1. Global Ocean Seasonal Data"Global Oceans"=select(global_summary, yr = year, temp_f = sst_f, anom_f) %>%mutate(rate = eq_global,rate_lab ="Global Ocean Warming"),# 2. Gulf of Maine Seasonal Data"Northeast US Shelf"=select(annual_summary, yr = year, temp_f = area_wtd_f, anom_f = area_wtd_anom_f) %>%mutate(rate = eq_all,rate_lab ="Northeast Shelf Warming")), .id ="region")# Compare warming ratesglobal_rate_plot <-ggplot(glob_plot_dat) +geom_line(data =filter(glob_plot_dat, region =="Northeast US Shelf"), aes(yr, anom_f), color ="gray20", linetype =3, linewidth =0.5, alpha =0.5) +geom_point(data =filter(glob_plot_dat, region =="Northeast US Shelf"), aes(yr, anom_f), color ="gray20", alpha =0.8) +# All year warming rate for GOMgeom_smooth(aes(yr, anom_f, color = rate), formula = y~x, linetype =1,method ="lm",se =FALSE) +# Annotation for Warming Rategeom_mark_ellipse(aes(filter = region =="Global Oceans"& yr ==2022,x = yr, y = anom_f,label = region), label.colour =gmri_cols("gmri blue"),con.colour =gmri_cols("gmri blue"),label.fill ="transparent", color ="transparent",label.fontsize =10) +# Annotation for Warming Rategeom_mark_ellipse(aes(filter = region =="Northeast US Shelf"& yr ==1999,x = yr, y = anom_f,label = region), label.colour =gmri_cols("orange"),con.colour =gmri_cols("orange"),label.fill ="transparent", color ="transparent",label.fontsize =10) +scale_color_gmri() +# Scalesscale_y_continuous(expand =expansion(add =c(1.5, 1.5)),labels =number_format(suffix =" \u00b0F")) +scale_x_continuous(limits =c(1981, 2025)) +labs(title =str_c("Northeast US Shelf & Global Ocean Warming Rates"),x ="Year",y ="Average Annual Temperature Anomalies",color ="1982-2022 Warming Rate") +theme(# legend.position = c(0.2, 0.85), legend.position ="bottom", legend.title =element_text(face ="bold"),legend.background =element_rect(fill ="transparent")) +guides(color =guide_legend(ncol =2, title.position ="top", title.hjust =0.5)) +coord_cartesian(clip ="off") +geom_logo(x_npc =0.9, y_npc =-0.15, logo_height =1, height_units ="cm") global_rate_plot```**Shift in Regional Ocean Circulation**Beginning around 2008, temperature in the region swung upward, and since 2010 most yearly temperatures have been over 1 degree above the hematological average of the previous 30-years. With some years experiencing temperatures more than 2 degrees C above that average.```{r}clim_avg <- annual_summary %>%filter(year %in%c(1982:2011)) %>%summarise(x =1982,xend =2011,avg_c =mean(area_wtd_sst, na.rm = T),avg_f =mean(area_wtd_f, na.rm = T),avg_anom_c =0,avg_anom_f =0) %>%mutate(across(where(is.numeric), round, 2))%>%mutate(across(where(is.numeric), round, 2))# average from 2010 onward:r2_avg <- annual_summary %>%filter(year %in%c(2010:2022)) %>%summarise(x =2010,xend =2022,avg_c =mean(area_wtd_sst, na.rm = T),avg_f =mean(area_wtd_f, na.rm = T),avg_anom_c =mean(area_wtd_anom, na.rm = T), avg_anom_f =mean(area_wtd_anom_f, na.rm = T)) %>%mutate(across(where(is.numeric), round, 2))#### Annual Trend Plot - regime shift ####ggplot(data = annual_summary, aes(year, area_wtd_anom)) +# Overlay yearly means and lines connecting themgeom_col(aes(fill ="Yearly Anomaly"), size =0.75) +scale_fill_manual(values =c("Yearly Anomaly"="gray90")) +# Add regression lines for shifted baselinesgeom_textsegment(data = clim_avg,aes(x = x, xend = xend, y = avg_anom_c, yend = avg_anom_c),size =4, linewidth =1,label ="1982-2011 Climatological Average",linetype =1, linewidth =1, color =gmri_cols("gmri blue"),hjust =0.50, vjust =-1.2) +geom_textsegment(data = r2_avg,aes(x = x, xend = xend, y = avg_anom_c, yend = avg_anom_c),size =4, linewidth =1,label ="2010-2022 Average",linetype =1, linewidth =1, color =gmri_cols("orange"),#boxlinetype = "dotted", boxlinewidth = 0.5,hjust =0.5, vjust =-1.2) +# Axesscale_y_continuous(sec.axis =sec_axis(trans =~as_fahrenheit(., data_type ="anomalies"),labels =number_format(suffix =" \u00b0F")),labels =number_format(suffix =" \u00b0C")) +# labels + themelabs(title ="Recent Temperatures Indicative of a Regime Shift",x ="", y ="Sea Surface Temperature Anomaly",caption =paste0("Anomalies calculated using 1982-2011 reference period.")) +# themetheme(legend.title =element_blank(),legend.position =c(0.85, 0.1),legend.background =element_rect(fill ="transparent"),legend.key =element_rect(fill ="transparent", color ="transparent")) ```The magnitude of this regime shift was `r r2_avg$avg_anom_c`C, or `r r2_avg$avg_anom_f`F. Bringing average yearly temperatures up from the 30-year average of `r clim_avg$avg_c`C to `r r2_avg$avg_c`C, or `r clim_avg$avg_f`F to `r r2_avg$avg_f`F.#### Seasonal Patterns```{r, fig.height=7}# Doing seasons by meteorological Definitionsquarter_summary <- region_timeseries %>%filter(year >=1982) %>%group_by(year = season_yr, season_eng) %>%summarise(sst =mean(sst, na.rm = T),sst_anom =mean(sst_anom, na.rm = T), area_wtd_sst =mean(area_wtd_sst, na.rm = T),area_wtd_anom =mean(area_wtd_anom, na.rm = T),.groups ="drop") %>%mutate(season_eng =factor(season_eng, levels =c("Winter", "Spring", "Summer", "Fall")))# Plotquarter_summary %>%ggplot(aes(year, area_wtd_anom)) +geom_col(fill ="gray60") +geom_smooth(method ="lm", aes(color ="Regional Trend"),formula = y ~ x, se = F, linetype =1) +stat_poly_eq(formula = y ~ x,color =gmri_cols("orange"),aes(label =paste(..eq.label..#, ..rr.label.., sep = "~~~" )),parse = T) +scale_color_manual(values =c("Regional Trend"=as.character(gmri_cols("orange")))) +labs(x ="", y ="Sea Surface Temperature Anomaly",caption ="Regression coefficients reflect annual change in sea surface temperature.") +scale_y_continuous(sec.axis =sec_axis(trans =~as_fahrenheit(., data_type ="anomalies"),labels =number_format(suffix =" \u00b0F")),labels =number_format(suffix =" \u00b0C")) +theme(legend.title =element_blank(),legend.position ="none",legend.background =element_rect(fill ="transparent"),legend.key =element_rect(fill ="transparent", color ="transparent")) +facet_wrap(~season_eng, ncol =1)```#### Overall Temperature Increase```{r, fig.height=2.25}dat_region <- annual_summary %>%filter(year %in%c(1982, 2022)) dat_global <- global_summary %>%filter(year %in%c(1982, 2022)) dat_list <-list(dat_region, dat_global) %>%setNames(c(tidy_name, "Global Oceans"))dat_combined <-bind_rows(dat_list, .id ="Area") %>%select(Area, area_wtd_sst, year) %>%pivot_wider(names_from = year, values_from = area_wtd_sst)ggplot(dat_combined, aes(x =`1982`, xend =`2022`, y =fct_rev(Area))) +geom_dumbbell(colour ="lightblue", colour_xend =gmri_cols("gmri blue"), size =3, dot_guide =TRUE, dot_guide_size =0.5) +geom_text_repel(aes(x =`2022`, y =fct_rev(Area), label =str_c("+", round(`2022`-`1982`, 2), " C")),color =gmri_cols("gmri blue"),vjust =4,hjust =0,segment.size =0.5,seed =123) +labs(x ="Sea Surface Temperature", title ="Change in Sea Surface Temeprature - 1982-2022", y ="") +scale_x_continuous(sec.axis =sec_axis(trans =~as_fahrenheit(., data_type ="temperature"),labels =number_format(suffix =" \u00b0F")),labels =number_format(suffix =" \u00b0C"))```### Marine HeatwavesA marine heatwave is defined as a situation when seawater temperatures exceeds a seasonally-varying threshold (usually the 90th percentile) for at least 5 consecutive days. Successive heatwaves with gaps of 2 days or less are considered part of the same event. The heatwave threshold used below was 90%. The heatwave history for `r tidy_name` is displayed below:**Note:** For the figures below heatwave events were determined using the methods of Hobday et al. 2016 and implemented using the R package [{heatwaveR}](https://robwschlegel.github.io/heatwaveR/).#### Heatwave Events```{r, fig.height = 6}# Plot heatmapheatwave_heatmap <-heatwave_heatmap_plot(hw_dat = grid_data, temp_units ="C")# Assemble piecesheatwave_heatmap```#### Heatwave Trends```{r, fig.width=8}#### Annual Heatwave Summary Details# number of heatwaves# average heatwave duration# remove NA as a distinct heatwave numberwave_summary <- grid_data %>%group_by(year(time), mhw_event_no) %>%summarise(total_days =sum(mhw_event, na.rm = T),avg_anom =mean(sst_anom, na.rm = T),peak_anom =max(sst_anom, na.rm = T),.groups ="drop") %>%drop_na() %>%group_by(`year(time)`) %>%summarise(num_waves =n_distinct(mhw_event_no),avg_length =mean(total_days, na.rm = T),avg_intensity =mean(avg_anom, na.rm = T),peak_intensity =max(peak_anom, na.rm = T),.groups ="drop") %>%rename(year =`year(time)`)#### Plotting# number of heatwaveshw_counts <-ggplot(wave_summary, aes(y = year, x = num_waves)) +geom_segment(aes(yend = year, xend =0), color =gmri_cols("gmri blue")) +geom_point(color =gmri_cols("gmri blue")) +labs(x ="Event Count", y ="")# average durationhw_lengths <-ggplot(wave_summary, aes(y = year, x = avg_length)) +geom_segment(aes(yend = year, xend =0), color =gmri_cols("orange")) +geom_point(color =gmri_cols("orange")) +labs(x ="Event Duration", y ="")# avg temphw_temps <-ggplot(wave_summary, aes(y = year, x = avg_intensity)) +geom_segment(aes(yend = year, xend =0), color =gmri_cols("green")) +geom_point(color =gmri_cols("green")) +labs(x ="Avg Anomaly", y ="")# peak temphw_peaks <-ggplot(wave_summary, aes(y = year, x = peak_intensity)) +geom_segment(aes(yend = year, xend =0), color =gmri_cols("teal")) +geom_point(color =gmri_cols("teal")) +labs(x ="Peak Anomaly", y ="")(hw_counts | hw_lengths | hw_temps | hw_peaks) +plot_annotation(title ="Heatwaves Frequent and Powerful in Recent Years")```# 2: `r current_yr` Observations### Temperatures of the Last Full Year (Interactive)Daily temperatures for the last 365 days are displayed below with reference to how they fall against the marine heatwave and cold spell thresholds.```{r}# # Option 1: Grab last 365 days# # # Grab data from the most recent year through present day to plot# last_year <- Sys.Date() - 365 #1 year from current date# last_yr_heatwaves <- region_hw %>% filter(time >= last_year)# last_yr <- year(last_yr)``````{r}# # Option 2: Grab last Full Year:# last_year <- Sys.Date() - 365 #1 year from current date# # # wind it back to first day of the year# last_year <- last_year - yday(last_year) + 1 # # # Filter out:# last_yr_heatwaves <- region_hw %>% filter(year(time) == year(last_year))# last_yr <- year(last_yr)``````{r}# Option 3: Last complete yearlast_year <-2022last_yr_heatwaves <- region_hw %>%filter(year(time) == last_year)``````{r}# Get number of heatwave events and total heatwave days for last year# How many heatwave events:# How many heatwave eventsnum_events <-max(last_yr_heatwaves$mhw_event_no, na.rm = T) -min(last_yr_heatwaves$mhw_event_no, na.rm = T) +1# Data for only the current yearthis_yr_hw <- region_hw %>%filter(year(time) ==2022)#this_yr_hw <- region_hw %>% filter(year(time) == year(Sys.Date()))# Number of heatwave eventsnum_hw_days <-sum(this_yr_hw$mhw_event, na.rm = T)``````{r mhw plotly}# Plot the interactive timeserieslast_yr_heatwaves %>%plotly_mhw_plots()```### `r year(Sys.Date())` TemperaturesTemperatures for the current year can be seen against the same thresholds with the following plot:```{r mhw ggplot}# Make Plothw_temp_p <-year_hw_temps_two(year_hw_dat = this_yr_hw, temp_units ="C")# Show Figurehw_temp_p```### `r year(Sys.Date())` AnomaliesChanging from absolute temperatures to anomalies reveals the degree to which this year is departure from what we would have expected from the climatology of 1982-2011's temperatures.```{r heatwave anomaly timeline}# Show figurehw_anom_p <-year_hw_anoms_two(year_hw_dat = this_yr_hw, temp_units ="C")hw_anom_p```### `r current_yr` in ContextWhen `r current_yr` is overlaid against the daily values and climatological average for the region, we can see how temperatures compare to the history of the region.```{r}# Plot this year among all previousthis_yr_hw <- this_yr_hw %>%mutate(yday =yday(time),flat_date =as.Date("2000-01-01") + yday -1)# Make a plotggplot(data = one_yr, aes(x = flat_date)) +geom_point(data = grid_data, aes(y = sst, color ="Daily Temperature"), alpha =0.6) +geom_textline(aes(y = seas, label ="Climatological Average", color ="Climatological Average"),linetype =1, size =4, hjust =0.5) +geom_line(data = this_yr_hw, aes(flat_date, sst, color ="2022 Temperatures"), size =1) +scale_x_date(date_labels ="%b", date_breaks ="1 month", expand =expansion(add =c(0,0))) +scale_color_manual(values =c("Daily Temperature"="gray80","Climatological Average"="gray10","2022 Temperatures"=gmri_cols("orange") ) ) +scale_y_continuous(labels =number_format(suffix =" \u00b0C")) +labs(x ="Day of the Year", y ="Sea Surface Temperature",title =str_c("2022 Temperatures in Context of the Last 40-Years"),subtitle ="Current year conditions suggest a persistance of warmer regime") +#theme_gmri() +theme(legend.position ="bottom",legend.key.height =unit(.5, "lines"),legend.key.width =unit(5, "lines"),plot.margin =unit(c(.5,1,.3,.5), "cm"),plot.title.position ="plot",legend.margin =margin(c(7,0,0,0)),legend.justification ="center") +guides(color =guide_legend(title ="",title.hjust =0.5,nrow =1,title.position ="left",label.position ="top",override.aes =list(shape =c(16, NA, NA),linetype =c(0, 1, 1),alpha =c(1, 1, 1),label =c("", "", "")),byrow = T))```# 3: Global Temperature MapsThe 2022 global sea surface temperature anomalies have been loaded and displayed below to visualize how different areas of the ocean experience swings in temperature.For perspective on where excess heat is distributed around the world, here first are maps of anomalies at a global perspective:```{r}# Access information to netcdf on boxnc_year <-"2022"anom_path <-str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", nc_year, ".nc")# Load 2022 as stackanoms_2022 <-stack(anom_path)```::: panelset::: panel## Annual Average {.panel-name}```{r}# Make Annual Averageann_anom_ras <-calc(anoms_2022, fun = mean, na.rm = T)# Color scale titlesst_lab <-"Sea Surface Temperature Anomaly"# Color limit for palettestemp_limits <-c(-2, 2)temp_breaks <-c(temp_limits[1], temp_limits[1]/2, 0, temp_limits[2]/2, temp_limits[2])temp_labels <-str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")# Make st objectann_anom_st <-st_as_stars(rotate(ann_anom_ras)) # Build Mapggplot() +geom_stars(data = ann_anom_st) +geom_sf(data = world, fill ="gray30", color ="white", size = .25) +scale_fill_distiller(palette ="RdBu", na.value ="transparent", limit = temp_limits, oob = scales::squish,breaks = temp_breaks, labels = temp_labels) +guides("fill"=guide_colorbar(title = sst_lab, title.position ="right", title.hjust =0.5,barheight =unit(3, "inches"), frame.colour ="black", ticks.colour ="black")) +coord_sf(expand = F) +map_theme() +theme(legend.position ="right", legend.title =element_text(angle =90)) +labs(title =str_c("Global Temperature Anomalies - ", nc_year))```:::::: panel## Winter {.panel-name}```{r}#### Making Averages for Seasons# Get previous year winterlast_nc_yr <-as.numeric(nc_year) -1last_yr_anom_path <-str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", last_nc_yr, ".nc")# Load 2020 as stack to get the full winterlast_yr_anoms <-stack(last_yr_anom_path)# Join to 2022anoms_double <-stack(list(last_yr_anoms, anoms_2022))# Drop december 2022 so its not influencing the Winter of 2021-2022dec_2022 <-"X2022.12"not_dec21 <-which(str_sub(names(anoms_double), 1, 8) != dec_2022)anoms_nodec <- anoms_double[[not_dec21]]# Set up list to use map()month_key <-list("Winter"=c("12", "01", "02"),"Spring"=c("03", "04", "05"),"Summer"=c("06", "07", "08"),"Fall"=c("09", "10", "11"))# Get mean anoms across seasonsseason_stacks <-map(month_key, function(mon){# Get names from the stack stack_names <-names(anoms_nodec) stack_months <-str_sub(stack_names, 7,8)# layers with correct months: in_season <-which(stack_months %in% mon)# Get mean across those months season_mean <-calc(anoms_nodec[[in_season]], mean, na.rm = T)# season_mean <- st_as_stars(rotate(season_mean))return(season_mean)})``````{r}# Plotting the seasonsseas_anom_maps <-function(season, lab_years, temp_lim =5){#temp_limits <- max(abs(values(season_stacks[[season]])), na.rm = T) * c(-1,1) # Center the color scale with Dynamic Limits temp_limits <-c(-1,1) * temp_lim temp_breaks <-c(temp_limits[1], temp_limits[1]/2, 0, temp_limits[2]/2, temp_limits[2]) temp_labels <-str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")# Make map seas_map <-ggplot() +geom_stars(data =st_as_stars(rotate(season_stacks[[season]]))) +geom_sf(data = world, fill ="gray30", color ="white", size = .25) +scale_fill_distiller(palette ="RdBu", na.value ="transparent", limit = temp_limits, oob = scales::squish,breaks = temp_breaks, labels = temp_labels) +guides("fill"=guide_colorbar(title = sst_lab, title.position ="right", title.hjust =0.5,barheight =unit(3, "inches"), frame.colour ="black", ticks.colour ="black")) +coord_sf(expand = F) +map_theme() +theme(legend.position ="right", legend.title =element_text(angle =90)) +labs(title =str_c(season, ": ", lab_years))return(seas_map)}``````{r}seas_anom_maps("Winter", "2021-2022", temp_lim =2)```:::::: panel## Spring {.panel-name}```{r}seas_anom_maps("Spring", "2022", temp_lim =2)```:::::: panel## Summer {.panel-name}```{r}seas_anom_maps("Summer", "2022", temp_lim =2)```:::::: panel## Fall {.panel-name}```{r}seas_anom_maps("Fall", "2022", temp_lim =2)```::::::# 4: Regional Temperature MapsThe regional views of 2022's sea surface temperature anomalies have been loaded and displayed below to visualize how localized differences changed throughout the year.```{r}# # Expand the area out to see the larger patternscrop_x <- crop_x +c(-2.5, 3.5)crop_y <- crop_y +c(-1.5, 1)# Make a new extent for croppingregion_extent_expanded <-st_sfc(st_polygon(list(rbind(c(crop_x[[1]], crop_y[[1]]), c(crop_x[[1]], crop_y[[2]]), c(crop_x[[2]], crop_y[[2]]), c(crop_x[[2]], crop_y[[1]]), c(crop_x[[1]], crop_y[[1]])))))# convert to sfregion_extent_expanded <-st_as_sf(region_extent_expanded)```::: panelset::: panel## Annual Average {.panel-name}```{r}# Mask the annual averagereg_ann_anom <-mask_nc(ras_obj = ann_anom_ras, mask_shape = region_extent_expanded)# Color limit for palettestemp_limits <-c(-5, 5)temp_breaks <-c(temp_limits[1], temp_limits[1]/2, 0, temp_limits[2]/2, temp_limits[2])temp_labels <-str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")# Build Mapggplot() +geom_stars(data =st_as_stars(reg_ann_anom)) +geom_sf(data = new_england, fill ="gray30", color ="white", size = .25) +geom_sf(data = canada, fill ="gray30", color ="white", size = .25) +geom_sf(data = greenland, fill ="gray30", color ="white", size = .25) +scale_fill_distiller(palette ="RdBu", na.value ="transparent", limit = temp_limits, oob = scales::squish,breaks = temp_breaks, labels = temp_labels) +guides("fill"=guide_colorbar(title = sst_lab, title.position ="right", title.hjust =0.5,barheight =unit(3, "inches"), frame.colour ="black", ticks.colour ="black")) +map_theme() +theme(legend.position ="right", legend.title =element_text(angle =90)) +coord_sf(xlim = crop_x, ylim = crop_y, expand = F) +labs(title =str_c("Regional Temperature Anomalies - ", nc_year))```:::::: panel## Winter {.panel-name}```{r}# Mask the Seasonsseasons_masked <-map(season_stacks, mask_nc, region_extent_expanded)# Plot the seasonsplot_masked_season <-function(season, year_lab, temp_lim =5){# Grab Season reg_seas_anom <- seasons_masked[[season]]# Get temp limits# temp_limits <- max(abs(values(reg_seas_anom)), na.rm = T) * c(-1,1) # Dynamic Limits# Color limit for palettes temp_limits <-c(-1, 1) * temp_lim temp_breaks <-c(temp_limits[1], temp_limits[1]/2, 0, temp_limits[2]/2, temp_limits[2]) temp_labels <-str_c(c(str_c("< ", temp_limits[1]), temp_limits[1]/2, 0, temp_limits[2]/2, str_c("> ", temp_limits[2])), "\u00b0C")# Build Map seas_map <-ggplot() +geom_stars(data =st_as_stars(reg_seas_anom)) +geom_sf(data = new_england, fill ="gray30", color ="white", size = .25) +geom_sf(data = canada, fill ="gray30", color ="white", size = .25) +geom_sf(data = greenland, fill ="gray30", color ="white", size = .25) +scale_fill_distiller(palette ="RdBu", na.value ="transparent", limit = temp_limits, oob = scales::squish,breaks = temp_breaks, labels = temp_labels) +guides("fill"=guide_colorbar(title = sst_lab, title.position ="right", title.hjust =0.5,barheight =unit(3, "inches"), frame.colour ="black", ticks.colour ="black")) +map_theme() +theme(legend.position ="right", legend.title =element_text(angle =90)) +coord_sf(xlim = crop_x, ylim = crop_y, expand = F) +labs(title =str_c(season, " - ", year_lab))return(seas_map)}``````{r}plot_masked_season("Winter", year_lab ="2021-2022", temp_lim =3.5)```:::::: panel## Spring {.panel-name}```{r}plot_masked_season("Spring", year_lab ="2022", temp_lim =3.5)```:::::: panel## Summer {.panel-name}```{r}plot_masked_season("Summer", year_lab ="2022", temp_lim =3.5)```:::::: panel## Fall {.panel-name}```{r}plot_masked_season("Fall", year_lab ="2022", temp_lim =3.5)```:::::: panel## Heatwave Progression {.panel-name}Looking specifically at the last heatwave event, we can step through how the event progressed over time, and developing pockets or warmer/colder water masses.```{r hw progression, eval = F}# Identify the last heatwave event that happenedlast_event <-max(region_hw$mhw_event_no, na.rm = T)# Pull the dates of the most recent heatwavelast_event_dates <- region_hw %>%filter(mhw_event_no == last_event) %>%pull(time)# Buffer the dates, start 7 days beforeevent_start <- (min(last_event_dates) -7)event_stop <-max(last_event_dates)date_seq <-seq.Date(from = event_start,to = event_stop,by =1)# Load the heatwave datesdata_window <-data.frame(time =c(min(date_seq) , max(date_seq) ),lon = crop_x,lat = crop_y)# Pull data off boxhw_stack <-oisst_window_load(data_window = data_window, anomalies = T, mac_os ="mojave")#drop any empty years that bug inhw_stack <- hw_stack[map(hw_stack, class) !="character"]##### Format the layers and loop through the maps ##### Grab only current year, format datesthis_yr <-stack(hw_stack)day_count <-length(names(this_yr))day_labs <-str_replace_all(names(this_yr), "[.]","-")day_labs <-str_replace_all(day_labs, "X", "")day_count <-c(1:day_count) %>%setNames(day_labs)# Progress through daily timeline to indicate heatwave status and severityhw_timeline <- region_hw %>%filter(time %in%as.Date(day_labs))``````{r hw progression animation, eval = F, animation.hook = 'gifski', fig.height=8, fig.width=9}#### Plot Settings:# Set palette limits to center it on 0 with scale_fill_distillerlimit <-c(max(values(this_yr), na.rm = T) *-1, max(values(this_yr), na.rm = T) )# Plot Heatwave 1 day at a time as a GIFday_plots <-imap(day_count, function(date_index, date_label) {# grab dates heatwaves_st <-st_as_stars(this_yr[[date_index]])#### 1. Map the Anomalies in Space day_plot <-ggplot() +geom_stars(data = heatwaves_st) +geom_sf(data = new_england, fill ="gray30", color ="white", size = .25) +geom_sf(data = canada, fill ="gray30", color ="white", size = .25) +geom_sf(data = greenland, fill ="gray30", color ="white", size = .25) +geom_sf(data =st_union(region_extent), color =gmri_cols("gmri blue"), linetype =2, size =1,fill ="transparent") +scale_fill_distiller(palette ="RdYlBu", na.value ="transparent", limit = limit) +map_theme(legend.position ="bottom") +coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +guides("fill"=guide_colorbar(title ="Sea Surface Temperature Anomaly \u00b0C",title.position ="top",title.hjust =0.5,barwidth =unit(3, "in"),frame.colour ="black",ticks.colour ="black")) # Set colors by name color_vals <-c("Sea Surface Temperature"="royalblue","Heatwave Event"="darkred","MHW Threshold"="coral3","Daily Climatology"="gray30")#### 2. Plot the day and the overall anomaly to track dates date_timeline <-ggplot(data = hw_timeline, aes(x = time)) +geom_line(aes(y = sst, color ="Sea Surface Temperature")) +geom_line(aes(y = hwe, color ="Heatwave Event")) +geom_line(aes(y = cse, color ="Cold Spell Event")) +geom_line(aes(y = mhw_thresh, color ="MHW Threshold"), lty =3, size = .5) +geom_line(aes(y = mcs_thresh, color ="MCS Threshold"), lty =3, size = .5) +geom_line(aes(y = seas, color ="Daily Climatology"), lty =2, size = .5) +scale_color_manual(values = color_vals) +# Animated Point / linegeom_point(data =filter(hw_timeline, time ==as.Date(date_label)),aes(time, sst, shape =factor(mhw_event)), color =gmri_cols("gmri blue"), size =3, show.legend =FALSE) +geom_vline(data =filter(hw_timeline, time ==as.Date(date_label)),aes(xintercept = time), color ="gray50", size =0.5,linetype =3,alpha =0.8) +labs(x ="", y ="",color ="",subtitle ="Regional Temperature \u00b0C",shape ="Heatwave Event") +theme(legend.position ="bottom")#### 3. Assemble plot(s) p_layout <-c(area(t =1, l =1, b =2, r =8),area(t =3, l =1, b =8, r =8))# plot_agg <- (date_timeline / day_plot) + plot_layout(heights = c(1, 3)) plot_agg <- date_timeline + day_plot +plot_layout(design = p_layout)return(plot_agg )})walk(day_plots, print)```::::::# 5: Ranking Warming RatesIf we look at the rates of change from 1982-2022 for each grid cell, rather than the observed temperature, it is possible to rank how hot each location on earth is warming relative to all the others.Once we have the rankings, we can then take the average ranking within the Gulf of Maine we can obtain the average warming rank for the area compared to the rest of the globe.```{r load warming rates}# 1. Warming Rates and Rankingsrates_path <-paste0(oisst_path, "warming_rates/annual_warming_rates")rates_stack_all <-stack(str_c(rates_path, "1982to2022.nc"), varname ="annual_warming_rate")ranks_stack_all <-stack(str_c(rates_path, "1982to2022.nc"), varname ="rate_percentile")``````{r mask warming rates}# Get the rank information that go with the original extent used by the timelinesranks_masked <-mask_nc(ranks_stack_all, region_extent)rates_masked <-mask_nc(rates_stack_all, region_extent)region_ranks <-get_masked_vals(ranks_masked, rates_masked)# Prep it for text input.avg_rank <- region_ranks$`Mean Rank`*100avg_rate <- region_ranks$`Mean Rate`low_rank <- region_ranks$`Min Rank`*100low_rate <- region_ranks$`Min Rate`top_rank <- region_ranks$`Max Rank`*100top_rank <-ifelse(top_rank ==100, "greater than or equal to 99.5", top_rank)top_rate <- region_ranks$`Max Rate````Based on data from 1982-2022, the warming rates of `r tidy_name` have been some of the highest in the world. The area as a whole has been increasing at a rate of `r avg_rate`$^{\circ}C/year$ which is faster than `r avg_rank`% of the world's oceans.Over that same period locations within the `r tidy_name` have been warming at rates as low as `r low_rate`$^{\circ}C/year$ and as rapidly as `r top_rate`$^{\circ}C/year$, corresponding to ranks as low as `r low_rank`% and as high as `r top_rank`%.Mapped below are the corresponding warming rates and their global rankings.```{r}# For the full map we mask again, but zoom out a little# Mask again using the expanded mask so we can zoom outranks_masked <-mask_nc(ranks_stack_all, region_extent_expanded)rates_masked <-mask_nc(rates_stack_all, region_extent_expanded)# Make stars objectrank_stars <-st_as_stars(ranks_masked)rates_stars <-st_as_stars(rates_masked)```### Warming Rates```{r}#### Rates Map ##### Make Contours if desired# rates_contour <- st_contour(x = rates_stars, na.rm = T, breaks = seq(0.85, 1, by = 0.02))# # rates map Original# rates_map <- ggplot() +# geom_stars(data = rates_stars) +# geom_sf(data = new_england, fill = "gray90", size = .25) +# geom_sf(data = canada, fill = "gray90", size = .25) +# geom_sf(data = greenland, fill = "gray90", size = .25) +# geom_sf(data = region_extent, # color = "black", # fill = "transparent", linetype = 2, size = 0.5) +# scale_fill_viridis_c(option = "plasma", na.value = "transparent") +# map_theme() +# coord_sf(xlim = crop_x, # ylim = crop_y, expand = F) +# guides("fill" = guide_colorbar(# title = "Annual Temperature Change \u00b0C/year",# title.position = "top", # title.hjust = 0.5,# barwidth = unit(2.5, "in"),# frame.colour = "black",# ticks.colour = "black"))# Look at the spread# hist(rates_stars$Annual.Sea.Surface.Temperature.Warming.Rate)# Reclassify to discrete binsrates_reclass <- rates_stars %>%mutate(warm_rate = Annual.Sea.Surface.Temperature.Warming.Rate,rates_reclass =case_when( warm_rate > .09~"> 0.1", warm_rate > .08~"0.08 - 0.1", warm_rate > .06~"0.06 - 0.08", warm_rate > .04~"0.05 - 0.06", warm_rate > .04~"0.04 - 0.05", warm_rate > .02~"0.02 - 0.04", warm_rate <0.2~"< 0.02",TRUE~"NA"),rates_reclass =factor(rates_reclass, levels =c("> 0.1", "0.08 - 0.1", "0.06 - 0.08", "0.05 - 0.06", "0.04 - 0.05", "0.02 - 0.04", "< 0.02","Outside Area Measured")) ) %>%select(rates_reclass)# Testing discrete scalesrates_map <-ggplot() +geom_stars(data = rates_reclass) +geom_sf(data = new_england, fill ="gray90", size = .25) +geom_sf(data = canada, fill ="gray90", size = .25) +geom_sf(data = greenland, fill ="gray90", size = .25) +geom_sf(data =st_union(region_extent), color ="black", fill ="transparent", linetype =2, size =0.5) +scale_fill_brewer(palette ="YlOrRd", na.value ="transparent", na.translate = F,direction =-1) +coord_sf(xlim = crop_x, ylim = crop_y, expand = F) +map_theme(legend.position ="bottom") +guides("fill"=guide_legend(title ="Annual Temperature Change \u00b0C/year",title.position ="top", title.hjust =0.5,label.position ="bottom", keywidth =unit(1.5, "cm"),reverse = T,nrow =1))rates_mapgrid::grid.raster(lab_logo, x =0.05, y =0.03, just =c('left', 'bottom'), width =unit(0.8, 'inches'))```### Warming Ranks```{r}#### Warming Ranks Map ##### Make Contours if desiredranks_contour <-st_contour(x = rank_stars, na.rm = T, breaks =seq(0.85, 1, by =0.02))# # Original Map# # ranks map# ranks_map <- ggplot() +# geom_stars(data = rank_stars) +# geom_sf(data = ranks_contour, fill = "transparent", color = "gray30", size = 0.1) +# geom_sf(data = new_england, fill = "gray90", size = .25) +# geom_sf(data = canada, fill = "gray90", size = .25) +# geom_sf(data = greenland, fill = "gray90", size = .25) +# geom_sf(data = region_extent, # color = "black", # fill = "transparent", linetype = 2, size = 0.5) +# scale_fill_viridis_c(option = "plasma",# na.value = "transparent",# limit = c(0.85, 1),# oob = scales::oob_squish) +# map_theme() +# coord_sf(xlim = crop_x, # ylim = crop_y, expand = F) +# guides("fill" = guide_colorbar(# title = "Global Percentile of Warming Rates",# title.position = "top", # title.hjust = 0.5,# barwidth = unit(2.5, "in"),# frame.colour = "black",# ticks.colour = "black")) +# labs(caption = "Ranking color scale truncated to display ranges of 0.85-1 # Lower values will display as 0.85 or 85%# Contour values every 2%")# Look at the spread# hist(rank_stars$Annual.Sea.Surface.Temperature.Warming.Rate.Rank)# reclassify to discrete binsranks_reclass <- rank_stars %>%mutate(warm_rank = Annual.Sea.Surface.Temperature.Warming.Rate.Rank,ranks_class =case_when( warm_rank > .999~"> 99.9%", warm_rank > .99~"99 - 99.9%", warm_rank > .98~"98 - 99%", warm_rank > .95~"95 - 98%", warm_rank > .90~"90 - 95%", warm_rank > .80~"80 - 90%", warm_rank < .80~"< 80%",TRUE~"NA"),ranks_class =factor(ranks_class, levels =c("> 99.9%", "99 - 99.9%", "98 - 99%", "95 - 98%", "90 - 95%", "80 - 90%", "< 80%")) ) %>%select(ranks_class)ranks_map <-ggplot() +geom_stars(data = ranks_reclass) +geom_sf(data = new_england, fill ="gray90", size = .25) +geom_sf(data = canada, fill ="gray90", size = .25) +geom_sf(data = greenland, fill ="gray90", size = .25) +geom_sf(data =st_union(region_extent), color ="black", fill ="transparent", linetype =2, size =0.5) +scale_fill_brewer(palette ="YlOrRd", na.value ="transparent", na.translate = F,direction =-1) +map_theme(legend.position ="bottom") +coord_sf(xlim = crop_x, ylim = crop_y, expand = F) +guides("fill"=guide_legend(title ="Global Percentile of Warming Rates",title.position ="top", title.hjust =0.5,label.position ="bottom", keywidth =unit(1.5, "cm"),reverse = T,nrow =1)) +labs(caption ="Pixels ranked globally based on annual warming rate.")# plot bothranks_mapgrid::grid.raster(lab_logo, x =0.05, y =0.03, just =c('left', 'bottom'), width =unit(0.8, 'inches'))```# 6: Shifting BaselinesIn 2022 NOAA is transitioning standard climatologies from the 30-year period of 1982-2011 to a new period spanning 1992-2020. Changes in climate regimes often does not result in a uniform upward or downward change that is consistent throughout the year.The plot below shows just how both the average temperature, as well as the annual highs and lows have shifted. When looking specifically at `r tidy_name` here is how the expected temperature for each day of the year has shifted.From this we can see that the Fall temperatures have risen more than those of the spring. There is also a large change in where the threshold for Marine Heatwave events sits, a consequence of exceptionally warm Fall temperatures becoming more common.```{r}# Run heatwave detection using new climate periodheatwaves_91 <-pull_heatwave_events(region_timeseries, threshold =90, clim_ref_period =c("1991-01-01", "2020-12-31")) %>%supplement_hw_data() %>%filter(doy !=366)# Subtract old from the newheatwaves_91 <- heatwaves_91 %>%mutate(clim_shift = seas - region_hw$seas,upper_shift = mhw_thresh - region_hw$mhw_thresh,lower_shift = mcs_thresh - region_hw$mcs_thresh)# Make arrows where we want to point at things:arrows <-tibble(x1 =as.Date(c("2000-07-15")),x2 =as.Date(c("2000-08-28")),y1 =c(1.25), y2 =c(1.15) )# Plot the differencesheatwaves_91 %>%filter(time >= last_year) %>%mutate(year =year(time),yday =yday(time),flat_date =as.Date(yday-1, origin = base_date)) %>%distinct(flat_date, .keep_all = T) %>%ggplot(aes(x = flat_date)) +geom_line(aes(y = clim_shift, color ="Mean Temperature Shift")) +geom_line(aes(y = upper_shift, color ="MHW Threshold Change")) +geom_line(aes(y = lower_shift, color ="MCS Threshold Change")) +geom_curve(data = arrows, aes(x = x1, y = y1, xend = x2, yend = y2),arrow =arrow(length =unit(0.08, "inch")), size =0.5,color ="gray20", curvature =-0.3) +annotate("text", x =as.Date("2000-06-01"), y =1.225, label ="Fall Extremes\nMore Frequent") +labs(x ="", y ="Shift in Expected Temperature \u00b0C",color ="") +theme(legend.position ="bottom") +scale_color_gmri() +scale_x_date(date_labels ="%b", date_breaks ="1 month", expand =c(0,0))````r insert_gmri_footer()`
